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Ultrafast diagnostics has boosted the importance of re- 
liable and efBcient simulations of quantum many-electron 
dynamics. The out-of-equilibrium evolution of systems 
ranging from small molecules to bulk materials has to be 
described on the fundamental scale over time intervals 
that exceed the relaxation time. The leading candidate 
for yielding tractable algorithms is the time-dependent 
generalization of the well-established Density Functional 
Theory (DFT) [l|: the Time-Dependent Density Func- 
tional Theory (TDDFT) [2-11]. To gain the status of a 
systematic and controlled description, however, TDDFT 
requires a firm footing. Its importance has been high- 
lighted by the number and diversity of TDDFT applica- 
tions 0, 01 and its potential impact has been heightened 
by the user-accessibility of numerical simulation tools ^ . 
But despite the continuous advances made since the sem- 
inal work of Runge and Gross Q, the foundations of 
TDDFT are still not as firmly established as those of 
DFT. This allowed some authors to openly cast doubt 
on the integrity of the approach 0-11] and to stand by 
their claim despite rebuttal from peers [lo| . A rigorous 
basis must be developed if TDDFT is to become an un- 
equivocal cornerstone of future quantum simulations. 

An important question that has remained open to this 
day is whether TDDFT can be rigorously formulated by 
means of a variational principle similar to that underly- 
ing DFT. Given the unifying and constructive role played 
by variational formulations in Physics and, especially, in 
DFT, one can expect that TDDFT would certainly ben- 
efit from such a formulation. In ground-state DFT 
the density no (r) of a system confined by a static poten- 
tial vo{r) can be obtained as the global minimum of the 
energy functional E[n] — F[n] + J drn{r)vo{r) by solving 



SE[n] 



Sn{r) 



0. 



where F[n\ is universal, i.e. independent of vq, 



F[n] = inf Tr 



sup 



E[v] 



D[T + W 

drn{r)v(r) 



(1) 

(2) 
(3) 



with T and W the kinetic and electron-electron energy 
operators. In Eq.([2]), the infinum is searched over all A^- 
particle density matrices D which yield the prescribed 
function n(r) via 



TrD7i(r) = n(r) 



(4) 



'Electronic address: ' daligaul®lanl.gav] 



This is known as the Levy-Lieb constrained-search pro- 
cedure, in honor to Levy and Lieb's contributions to sub- 
stantiate the original proposal of Hohenberg and Kohn 
[l2,[I3l- Equivalently, F[n] can be obtained as the Leg- 
endre transform ^ of the ground-state energy E[v] = 
inf>p(^']i/ Ji;]]^') when the latter is considered as a func- 
tional of the external potential v{r). 

The establishment of a variational principle capable 
to underlie TDDFT has remained a challenging problem 
for many years that non only gives rise to some deep 
problem of analysis akin to those encountered in ground- 
state DFT but also to difficulties specific to the time- 
dependent case d, IT3 - l2Ci| . It was realized that the ini- 
tial proposal of Runge and Gross ^ based on the Dirac- 
Frenkel action principle is inadequate since it is incon- 
sistent with causality requirements. The problem arises 
since the Dirac-Frenkel principle wrongly fixes boundary 
conditions at both end points of the action functional, 
which is incompatible with the first-order character of 
the Schrodinger equation [2l|. The difficulty is not un- 
like that encountered in quantum many-body physics and 
field theories to extend ground-state perturbation theory 
to non-equilibrium states. There, a general and widely 
used solution was developed by Schwinger and Keldysh, 
and consists in extending the dynamics along a closed 
contour in time[22[. Using this artifice, van Leeuwen 
[m (and also jl6l]) suggested a route to define an varia- 
tional action principle that respects causality, but whose 
steps still remain to be justified. 

In this letter, we present a rigorous constrained-search 
variational formulation of TDDFT that extends the 
Levy-Lieb formulation of ground-state DFT to the time- 
dependent case. The TD particle density can be obtained 
as the stationary point of an effective action functional, 
which is itself defined from a constrained variation of the 
quantum action of Balian and Veneroni f^s'] (in place of 
the Frenkel-Dirac action); the approach builds on ideas 
originally developed in field theory [H, [2^ . In the spirit 
of Levy-Lieb, a universal action functional is defined and 
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is related to the functional introduced by van Leeuwen 
[l8| . The formulation is then used to justify the TD 
Kohn-Sham scheme. For clarity, we successively define 
our notations, state the main results and derive them. 

Definitions. We are interested in the dynamics over 
the time interval to < t < t f oi a. system of N electrons 
under the influence of an external, time-dependent scalar 
potential wo(r, t); we assume that the system evolves from 
a fixed, known initial state at to described by the density 
operator Dq; the particle density is denoted by no(r, t). 

To circumvent difficulties with causality, we regard the 
dynamics under investigation as a special case of a fam- 
ily of dynamics (i) along a closed-time contour ^(t) on 
C = [Ti,Tf] that monotonically goes from to to tj and 
monotonically back to to as r runs from to r/ , and (ii) 
governed by the Liouville equation. 



^^dlXrl ^ ^ H^{t),D{t] 
dr 



with b{n)=bo,{b) 



where Hu{t) = T + W + f drn{r)u{r , t) and u{r,T) is an 
external potential on C. With u = Mo(r,r) =vo{r,t{T)), 
u is equal on both branches of C, and the dynamics 
governed by Eq.® follows that of the system under in- 
vestigation along the forward branch and retraces back- 
ward the same trajectory along the return branch; there- 
fore, the particle density is also equal on both branches 
to the interacting density no(r, ^(t)). In general, how- 
ever, different potentials may act on each branch and, 
in turn, the forward and backward evolutions differ. 
For later reference, we define the evolution operator 

Uu{t) = Tks. -^^i ^^Yieve Tk is time-ordered 

chronological operator along the closed-time contour; we 
also introduce the Heisenberg representation of the par- 
ticle density along C, h^{Y,T) = ?7^7-^(r)n(r)J7t,(r); thus 
no(r,t(T)) = TrZ)onMo(r,T)- 

The Liouville equation ([5]) can be derived b y ex trem- 
izing the extended Balian-Veneroni functional [23| . 



ihiTrOrfVrj, - 1 



^ dO Y - - 1 \ / . , 



(6a) 
(6b) 



under arbitrary variations of the operators (i.e. matrix 
elements) Vr and Or, and subject to the boundary con- 
ditions 



Do 



Orf = 1 (identity operator). 



(7a) 
(7b) 



For our purposes, we have extended the original Balian- 
Veneroni action functional along the closed-time contour 
and added an insignificant constant ih to set its station- 
ary value to zero. Demanding Aui'D, O) to be stationary 
against arbitrary variations of Vr and Or not only gen- 
erates Eq.([5]) for D{t) but also a Liouville equation for 
0{t) with boundary conditions (j7bp at t/. However, the 
latter selects the constant solution 0{t) — 1 for all r £ C 
and one is left with Eq.([5]) only. Nevertheless the ap- 
parently superflous variational parameter Ot in Eq.Q is 
necessary to ensure that only one boundary condition be 
associated with the Liouville equation ([S]), in conformity 
with its first-order character and thus with causality re- 
quirements. As for the closed-time contour, it ensures 
causality of the effective action r[n] defined below. 

Main results. We separate them into four propositions. 
For ease of comparison. Table |T] gives the correspon- 
dence between the Levy-Lieb formulation of DFT and the 
present extension to TDDFT. Below, Jdr = J drt'ir). 
1- Constrained- Search Variational Principle: Consider 



DFT 


E[n] 


E[v] 


F[n] 


Eq.© 


Eq.® 


Eq.Q 


TDDFT 


r[n] 


A[v] 


F[n] 


Eq.(fTT1) 


Eq.(fT3l) 


Eq.© 



TABLE L Correspondence between DFT and TDDFT. 



the ejjfective action functional 



r[n] = s.v. AuoiV.O) 



(8) 



defined for all n e A/" (defined below) as the stationary 
value ("s.v.") of Aug when Vr and Or are varied subject 
(i) to the constraint of density. 



1 



[Vr.Or} n(r) 



n(r,T) , 



(9) 



(which upon integration implies TiVrOr = 1 ), (ii) to 
the contraint of phase. 



TrVr = 1 



TrC^ = 1 



(10) 



and (iii) to the boundary conditions ([7]). A/" is the 
set of generalized particle densities consisting of func- 
tions n{r,T) defined on R'^ x C that can be expressed 
as in Eq.® and such that ?T,(r, r^) = no(r,to) and 
J drn{r, r) = N. Then, the TD density no of the sys- 
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tern under investigation is a stationary point of r[n], 



ST[n] 



6n(r, t) 



0. 



(11) 



The formulation involves a two-stage, constrained search 
akin to the Levy-Lieb procedure: first an effective action 
r[n] is constructed from a constrained variation of an 
action, which is then stationarized to obtain no. How- 
ever, while in DFT the stationary value corresponds to a 
global minimum, it is difficult to characterize further the 
stationary point obtained using T[n] [IB]- That should 
not be a worry since problems of motion are generally 
not influenced by the specific extremum conditions; for 
instance, in classical mechanics, the least action principle 
is a stationary principle and the additional criteria of a 
true extremum are of interest only if stability is involved. 
2- Dual Representation and Universal Effective Action 
Functional. The effective action satisfies 



r[n] 



F[n] - dr J 



druo{r,t{T))n{r,T) (12) 



in terms of the universal (independent of uq) functional 



F[n] = s.v. 



+ / dr dru(r,T)n{r,T) 



(13) 



where A[u] = TrDoTxe R -^^^ dTH^(r) -g ^j^^ Leeuwen 
functional introduced in [isj (here extended to mixed 

states since [ll] assumes Dq = |^'o)(^'o|). 
3- Kohn-Sham approach. The functional F can be de- 
composed as F[n] = FKs[n] + Fex[n] such that (i) Fks 
is independent of the interparticle interactions, (ii) the 
electron-electron interaction contributions are contained 
in the excess term F^x, and (iii) the potential defined by 



UKS 



5FKs\n] 



Sn 



5F,. 



Uq 



Sn 



(14) 



is physical, i.e. equal on both branches, u/f5(r,T) 
VKs{'r,t{T)). The quantity TksM = FksM - 
J^' dr J drvKs{^,t{T))n{r,T) is the effective action of a 
noninteracting system of particles in the external poten- 
tial vks with the same density uq as the fully interacting 
system. Its dynamics is governed by the single-particle 
Liouville equation 

- [pV2m + VKsir, t), pKs] , (15) 
so that no{r,t) — TrpKs{t)n{r). 

4- Static limit: The present formulation can be applied 
to static external potentials, vo{r,t) = t'o(r), and equi- 
librium initial state Dq = e-''(*+^+^ at (in- 
verse) temperature /?. It provides an alternative to the 
Levy-Lieb formulation of DFT. The situation is not un- 
like the duality that occurs in the theory of equilibrium 
Green's functions between the Keldysh and the finite- 
temperature (Matsubara) approaches 0, [l^l . 



Proofs. 1- The proof proceeds by carrying out the 
s.v. calculation ([5]) explicitly using a Lagrange multiplier 
J(r, r) to enforce under independent variation of Vj- 
and Or, two equations emerge. 



ih 



1 dD{T) 

dr 

1 dd{T) 
V(t) dr 



K{r),D{T)\ - J drJ(r,r)i{^(r),n(r)}^ 

= with D{n) = A (16a) 

ff„„(T),0(r)] + 1 drJ(r,r)l{0(r),n(r)}^ 

= withO(T/) = /. (16b) 



These can be regarded as constrained Liouville equa- 
tions, where the free dynamics is subjected to an addi- 
tional non-mechanical "force" that maintains the aux- 
iliary condition ([9]) at all times along C. In what fol- 
lows, we first derive a set of relations implied by ()16|) 
and then apply them. As we shall see, at that stage, 
the Lagrange parameter is defined up to a purely time- 
dependent constant; to keep track of it, we replace J(r, r) 
by J(r,T) +w{t)/N in Eq. ([T6| . Then as may be verified 
be direct substitution [28|, the solution of Ea. ([TE)) is 



Dir) = {t)C{t, n)DoC{r, n)Uu,}{r)e'^ ''^^'\n) 
d{T) = Uu„{r)C{Tf,r)C{Tf,T)U-^\r)e-i^^''''-('^' 



f f d-rdrJ{r,T)n^g{r,T) 

with C(r, r') Tkc ^' ; we note 

that D{t) and 0{t) carry a phase related to w{t). Using 
Eqs.ldH) in TrP(T)d(T) ^ 1 yields 



dTw{T) = ihlnTrDaUi{n,Tf) = W[J] (18) 



with Ui{t,t') = Tife~^J'x''*^/'*"-"'('-'^)""o("-'^). Moreover, 
substituting Egs. pTl) in the density constaint © gives 



_ n(r,T): 



TiDoUiiTf^n) 



Tr£)oTif;7/(T/,T)n„o(r, T)Ui{T,Ti) , 



(19) 



which is independent of w{t). Incidentally, the right- 
hand side of Eq. ([T9l) is also the first derivative of W[J] 
defined by Eq. p^ . and therefore 



n(r, r) 



SW[J] 



(20) 



<5J(r,T) ■ 

Finally, using Eas. (|16p . the s.v. ([SJ can then be expressed 



r[7i] = -W[J] + J ' dr J dr J(r, T)n(r, 



(21) 



Using Eg. p^ . we find that for all purely TD quantities 
w{t), 

-W[J+w] + [{J + w)n^ ^W[J\ + [jn^ T[n] . 
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and therefore T[n] is well defined: it is single- valued de- 
spite the degeneracy in the Lagrange multiplier for a 
given n. The normalization constraint (jlOp allows one 
to fix w{t) in Eqs.([T71) to 



w{t) 



TtDoTk J drJ(r,T)n„„(r,r)e 



-iff drdrj (r,T)n„Q (r^r) 



TiDoTKe 



/ /drdr,/(r,r)ri„o(r,r) 



and therefore to define a one-to-one mapping between n 
and J (the J n mapping is defined through Eq. ([20| ). 
In particular, with n = rip, this gauge selects the natural 
solution J = and w = 0. Under these conditions, 
Eq. ([20l) says that n and J are conjugate variables and 
the relation (|2T|) for r[n] can be regarded as the Legendre 
transform of M^[J]; it follows that 



Sn{r, t) 



J{v,t). 



(22) 



The removal of constraints, which corresponds to the 
physical situation under consideration, is equivalent to 
J = 0; hence, from Eg. (|22|) . T[n] is stationary with re- 
spect to the variation of n around no. 

2- Using Uuo+j{Tf,Ti) = UuoiTf,Ti)Ui{Ti,Tf) and 
Uuo{Tf,Ti) = 1, we find 

W[J] = A[u =ua + J]. 
Using Eq. (PT|) . it follows (dropping variables) 



r[n] 



-iK] + I nJ 















[ nu* 










= s.v. 


-A[u] ~ 






U 







UqJ 



UqJ 



F[n] - J UqJ , 



which proves proposition 2. 

3- A complete proof is beyond the scope of the paper, 
and we hope to publish it elsewhere. We outline the main 
elements and point to relevant references for details. By 
treating the strength of the electronic Coulomb repulsion 
g = as an expansion parameter, the effective action 
F[n\ can be expressed as 



(23) 



k=0 



Different techniques can be used to obtain the coupling- 
constant expansion (^5)) from Eq. ([T^. such as the 



auxiliary field method [l6|, |29j or the inversion method 
The zeroth-order contribution F^s = 
F^'^\ which is the only remaining term in the limit 
5 — 0, corresponds to the universal effective action of 
a non-interacting system, while the higher-order correc- 
tions collected in Fex = X^fc^i ff*^^'''^'' ["■] contain the ef- 
fect of electron-electron interactions, including the TD 
Hartree and KS exchange contributions at lowest order. 
Then, Eq. (fT4|) is obtained by combining 6T/Sn\no = 
6F/6n\nf, — Uq — with F = F^s + Pex- From the sym- 
metry property of the effective action under exchanging 
the forward and backward components of n along C, the 
solution of SFKs/Sn evaluated at the physical density uq 
can be shown to be equal on both branches. Finally, from 
propositions 1 and 2, ri<-5[n] = FKsin] — J vks^ can be 
interpreted as the effective action s.v. Auj^g {T>, O). At 

n — hq, the stationary point implies the free Liouville 
equation, which, for a system of independent particles, is 
equivalent to the single-particle Liouville equation ([15]) 
for the one-particle density operator pks- For a pure 
initial state, the latter can be diagonalized and one ob- 
tains the more usual single-particle KS equations for the 
KS orbitals. We remark that, in the absence of a varia- 
tional principle, the TDKS scheme was either justified by 
assuming the so-called non-interacting w-representability 
of interacting densities, or searched by direct construc- 
tion of the KS system through solution of a non-trivial 
partial-differential equation for vks [H, • 

In conclusion we have described a constrained-search 
variational formulation of TDDFT. The key results are 
collected in the four propositions listed above and the 
main ingredients are conveniently compared with those 
of the Levy-Lieb constrained-search formulation in Table 
(U]). The approach integrates and extends ideas previously 
developed in other fields [H, [H, [H, (2^. A profound 
mathematical analysis of its underpinnings, well beyond 
the scope of the present paper, would be much desirable; 
we hope that the incompleteness of the results presented 
here will encourage others to pursue some of the ques- 
tions raised by them. We anticipate that many of the 
techniques developed in field theory, e.g. functional inte- 
gration, loop expansion, auxiliary field..., could be much 
valuable to develop controlled approximations based on 
r[n] and to provide additional insights on the TD Kohn- 
Shani potential. Finally, the approach can be expanded 
to TD-current-DF T p and can also be adapted to purely 
classical systems fM] by considering the classical coun- 
terparts of the different quantities (e.g., density opera- 
tors into distribution functions, commutator into Poisson 
bracket...) 

This work was carried out under the auspices of the 
National Nuclear Security Administration of the U.S. De- 
partment of Energy at Los Alamos National Laboratory 
under Contract No. DE-AC52-06NA25396. 
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